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Abstract 

In this paper we consider the Sturm-Liouville equation —y" + qy = Xy on the half line (0, oo) under 
the assumptions that x — is a regular singular point and nonoscillatory for all real A, and that 
cither (i) q is L\ near x = oo, or (ii) q' is L\ near oo with q{x) — > as x — > oo, so that there is 
absolutely continuous spectrum in (0,oo). Characterizations of the spectral density function for this 
doubly singular problem, similar to those obtained in |12j and |13j (when the left endpoint is regular) 
are established; corresponding approximants from the two algorithms in [T^] and [T3] are then utilized, 
along with Frobenius recurrence relations and piecewise trigonometric - hyperbolic splines, to generate 
numerical approximations to the spectral density function associated with the doubly singular problem 
on (0, oo). In the case of the radial part of the separated hydrogen atom problem, the new algorithms are 
capable of achieving near machine precision accuracy over the range of A from 0.1 to 10000, accuracies 
which could not be achieved using the SLEDGE software package. 
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1 Introduction 



In this paper we consider the Sturm-Liouville equation, 

- y" + q(x)y = Xy, (1.1) 

on (0, oo), with two singular endpoints on the half line (0, oo) under the same assumptions as in [101 114) 
under which x = is a regular singular point of type LC/N or LP/N (limit circle or limit point 
and nonoscillatory at x = for all real A) and x — oo is of type LP/O-N with cutoff A = (limit 
point and nonoscillatory at x = oo for A e (— oo,0) and oscillatory for A € (0, oo)); see [HI p. 114] 
for these definitions of endpoint classifications developed in connection with the SLEDGE software 
package. Under these assumptions, the spectrum is simple and the eigenfunction expansion associated 
with equation (|1.1[) has the general form 



/(*) = 



£ 

A„<0 



T(A) • <p{x,\)dp{\) 

f(t)4>(t,\n)dt 



</>(x, A„) 



T(X)-cf>(x,X) dp (A) 



(1.2) 



where 



T(X) := lim / f(x)(f>(x,X)dx, 
b ^°° Jo 



(1.3) 



and the solution </>(•, A) is a suitably normalized Frobenius solution near the regular singular point x = 0. 

If we make, in addition to the above assumptions, the more stringent assumptions posed in [12j (see 
Assumption 3 below), then there is absolutely continuous (a.c.) spectrum in (0,oo), and the spectral 
function p(X) in (|1.2I) is absolutely continuous on all closed intervals in (0, oo). 

The purpose of this paper is (i) to extend the analysis in [121 113] under suitable assumptions on 
q(x) to show that the spectral density function associated with (|1.2|) over the a.c. range (0, oo), that 
is, 



f(X)=p'(X), p(X) 



f(p) dp 



(1.4) 



can be represented for all A £ (0, oo) as (see Theorem 3 below) 



/(A) 



n[P(x, X)<j){x, A) 2 + Q(x, A)0(x, X)c/)'{x, X) + R(x, X)<f)'(x, A) 2 



(1.5) 



where (P(-, A), Q(-, A), R{-, X)) T is the unique solution of the initial value problem at x = oo 



dU 
dx 

lim 



d 

dx 



P 

Q 

R 



P(x,X) 
Q(x,X) 
R(x,X) 






' p ' 




Q 

R 



(1.6) 



Ae(0,oo), (1.7) 
and (ii) to extend the numerical algorithms from [T21[r3] for the computation of the spectral density 
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function /(A). We illustrate the new numerical algorithms on several examples, including the radial 
part of the separated hydrogen atom. For the first objective (i) we make use of the fact that /(A) 
is characterized as the boundary value of a suitable Titchmarsh-Weyl m-function by the Titchmarsh- 
Kodaira formula, which was recently established for such doubly singular problems in |16[ 114] (see 
Theorem 2 below). For the second objective (ii) we make use of exact Frobenius power series to 
estimate the solution 4>{x, A) and its derivative near x = 0, and then apply initial conditions at a 
suitable point x Q (X) > using values of (p(xo(X), A) and (f>'(xo(X), A) which can generally be computed 
to machine precision, so that numerical algorithms from |121 1 1 3) with the left endpoint regular can be 
adapted to approximate the right hand side of (ll.5[) ; this is done by shooting with piecewise trigonometric 
/ hyperbolic splines to compute the solution <f> of fll.10 and the solution (P, Q, R) T of (]1.6p at a suitable 
'matching' point x £ (xo(A),oo). 

The organization of topics needed in this paper to accomplish the above two objectives is as follows: 
In section 2 we give the main assumptions near x — from [10] fl4] , and the general forms of two linearly 
independent Frobenius power series solutions in all the cases we consider in this paper. In section 3 we 
list (without proof) the elementary results which relate and interconnect solutions of the Sturm-Liouville 
equation (jl.ll) with solutions of Appell's first order system f| 1 : none of these elementary results require 
any special assumptions on the potential q. In section 4 we add the main assumptions from [12 under 
which the initial value problem (ll.6[) - (|1.7[) has a unique solution in (0, oo), and reformulate in terms 
of solutions of Appell's system (|1.6I) results obtained by D.B. Pearson and his student Al-Naggar in 
[TJ [5] . This yields the spectral density function characterization (|1.5p in the relatively simple case of 
a regular left endpoint. In contrast to [TJ [2] we do not focus the analysis on the third order ordinary 
differential equation (see f|4. 13[) below) which is satisfied by the third component, R(x,X), of (|1.6p . but 
make use instead of many of the elegant formulas from section 2. In section 5 we generalize the methods 
of [U[2] to the doubly singular problem on (0, oo), when the m-function is defined relative to the suitably 
normalized Frobenius fundamental system as in [10j[l4]. This yields Theorem 3 below (a new result) 
in which the new Titchmarsh-Kodaira formula (see Theorem 2 below) for the (single) spectral density 
function associated with the doubly singular problem gets converted to the form (|1.5[) . In section 6 we 
list four test examples of equations on (0, oo) from [10] 114] for which explicit closed form formulas for 
the spectral density function were obtained. In section 7 we describe how the two different types of 
numerical algorithms from [12] and [13] (for cases involving a regular left endpoint) can be adapted, 
using appropriate heuristics, to yield new algorithms for computing the spectral density function when 
both endpoints are singular, which is done by utilizing the characterization (|1.5[) in an appropriate way. 
In section 8 we give numerical output showing that the new algorithms for doubly singular problems 
can achieve very high accuracy on the four test examples in section 6; we also give numerical ouput 
demonstrating convergence of our numerical approximations for a potential on (0, oo) from quantum 
chemistry where the potential has an infinite series representation satisfying all our assumptions. In 
section 9 we make use of our new code, AutoB, for the spectral density computation in order to generate, 
by quadratures, approximations to the spectral functions for the four examples in section 6, and give 
comparisons on timing and accuracy with the corresponding SLEDGE runs. Our main conclusion is 
that the new algorithms for doubly singular problems are very much superior to the older algorithms 
for doubly singular problems which were implemented in the SLEDGE software package. 
Remark. In our previous papers [T2J H5] the system (]1.6p was referred to as the "PQR equations" (our 
notation); however, the analysis leading to them (particularly the motivating property (]3.14l) ) was 
discovered by M. Appell [3] in 1880. Accordingly, we will henceforth refer to this first order system as 
the Appell equations. 

2 Suitably Normalized Frobenius Solutions 

In this section we repeat the basic definitions and some of the elementary properties of the suitably 
normalized Frobenius solutions which were introduced in Fulton |10] and Fulton and Langer [14] . It 
will be noted that, in most of the common cases (Bessel quations, Confluent Hypergeometric equations, 
Whittaker equations) rather standard normalizations of well known special functions have to be aban- 
doned in order to achieve the desired analytic properties of the Frobenius solutions (particularly, entire 
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behaviour in A) needed for the fundmental definition of a single Titchmarsh-Weyl m-function, the cor- 
responding scalar spectral function, and for determination of eigenfunction expansions of the problems 
considered in this paper (all of which have simple spectrum) . 

We consider in this paper the Sturm-Liouville equation (jl.lj) on the half line (0, oo) under the following 
assumptions (see [TU1 H^] ): 

Assumption 1: Near x = : 

Case I: For all x G (0,oo), 

oo 

q(x) = + — + q n+2 x n , q n real for all n, (2.1) 
x z x ■* — ' 

n— 

where the series is convergent in (0, oo), and where 

— j < <7o < oo. an d qo,qi not both zero. (2-2) 

or 

Case II: There exists a > such that q(x) is given by (|2.1[) for x E (0, a] where the series is convergent 
in (0, a] and where (|2.2[) holds, and in the interval [a, oo) we have q G L^Ja, oo), 

and 

Assumption 2: Near x = oo : 



lim q{x) = 0, (2.3) 

a;— >oo 

The assumptions near x = ensure that the indicial roots near the regular singular point x = are both 
real; it follows that the endpoint x = is either LC/N or LP/N, and the assumption near x = oo 
ensure that the endpoint x = oo is LP/O-N with cutoff A = in the terminology of [15]. Under the 
above assumptions it was proved in [101 Theorems 4.2,4.3,5.3,5.4] and [HJ Theorem 4.5] that the 
eigenfunction expansion associated with (|l.lj) (in both the LP and LC cases at x = assumes the 
form (|1.2j) . with a suitably normalized Frobenius solution </>(■, A). 

We now give the formulas for all cases of Frobenius solutions which can occur at x — under the 
assumptions (|2 . 1 [) - (|2 . 2[) ; these are the solutions which were utilized in [TUltH]. The indicial equation for 
the R.S.P. x = for the Sturm-Liouville equation (|l.ip with potential (|2.1I) . 

- V"{x) + [ % + ?L + £ qn+2X ™ ) y ( x ) = X y(x), x G (0, oo). (2.4) 



is 



X x 

n=0 



..2 .; I ,;2 1 » ' / 1 \ \ / /l 



H-r-gb=T--r-(^--l = fr-(- + i/ll.(r-l--t/J] (2.5) 



where we have set 



q Q = v 2 - ^v>Q 



for convenience. This gives rise to the following cases of Frobenius solutions: 

Case I: -\ < q < oo, q Q = v 2 - \ ^ M ^ -1 , M = 1,2,... .(This is Case I A in [TP]). In this case, 

Ittfo A) = ^1 + £ a n (A)a; n ^ , (2.6) 
y 2 (x, A) = a;'-" ^1 + 6„(A)x"^ , (2.7) 
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where a n (X), b n (X) are polynomials in A of degree [5] , and 

W x ( yi (-,X),y 2 (-,X)) = -2v. (2.8) 

Case II A: q = M odd: M = 2£+ 1, £ = 0,1,..., that is, q =£(£ + 1). 

In [TU] this is Case IC for M odd and it includes Case II (for I = 0). In this case, 

yi(x,X) = x e+1 ^ + £ a n (X)x n ^j , (2.9) 

y 2 (x,X) =K e (\)y 1 (x;\)\nx + x- t f 1 + £ <UA)z n j , (2.10) 

where a„(A), d n (X) are polynomials in A of degree , iff (A) is a polynomial of degree ^, and 

W x {yi(;X),y 2 (-,X)) = -(2£+l). (2.11) 
Case II B: q = ^f 1 , M even: M = 2N, N = 0, 1, . . . , that is, q = N 2 - \. 
In [10] this is Case IC for M even and it includes Case IB (for N = 0). In this case, 

yi (x,X)=x?+ N ^l + £a„(A)x^ , (2.12) 

OO 

V2(x,\) = yi (x,X)lnx + J2 d n (X)x? +n , if N = 

n=l 

y 2 (x,A) = ^Ar(A)yi(x;A) lux + a:*"* ^ + £ ^n(A)a;"J , if N > 1 (2.13) 

where a„(A), d„(A) are polynomials in A of degree [^] , K N (X) is a polynomial of degree N, and 

f -2iV if N>1, 
W x ( Vl (;X ,y 2 (-,A -i ., ," n (2-14) 

I 1 if N = 0. 

In each of the above cases the first Frobenius solution yi(-,X) is the principal solution at x = for 
all A e ( — oo, oo). The Frobenius solutions as normalized above satisfy the following properties (see [TU1 
Theorem 2.1]): 

(i) yi(x, •), 2/2(2;, •) and their derivatives are entire functions for each x € (0, 00) and satisfy for all 
A £ G, x £ (0, 00) the relations 



Vi (x, A) = yi(x, A), y- (z, A) = y^(x, A), i=l,2. 
(ii) yi(-, A) £ £2(0, xo) for < xo < 00 and for all A € (E. 
(iii) W x (yx(-,X), y 2 (-, A)) = C ^ where C £ R, independent of A. 

It follows that a fundamental system of solutions near x = which is entire in A and satisfies property 
(i), together with the normalization 

W x (ct>{;X),6(;X)) = 1 for all X £ (D, (2.15) 

can be selected by taking 

<f>{x, A) := yi (x, A), 9{x, A) := y 2 (x, X)/C, (2.16) 

where C is the real constant in (|2.8I) . (|2.1ip . or (|2.14p . In the LP cases at x = only </>(•, A) 
satsifies the property (ii) of square integrability near x = 0, so in all LP cases it is the first Frobenius 
solution which is used to write the eigenfunction expansion in the form (11.2[) . 

The LC cases at x = are Case I with q Q £ (-±, 0) U (0, §), Case IIA with q = 0(£ = 0), and 
Case IIB with q = \(N = 0). In this paper we limit our consideration of LC boundary conditions at 
x = to the Friedrichs LC boundary condition (sec (|5.2p below); in this case it is the first Frobenius 
solution cj> (the principal solution) which is selected and used in the eigenfunction expansion (|1.2[) . 



5 



3 Preliminaries 



In this section we collect together some useful results which relate solutions of the Sturm-Liouville equa- 
tion to solutions of the companion first order system (|1.6p . The proofs of these results (though 
sometimes tedious) require only straightforward algebraic manipulation making use of these two equa- 
tions, and no special assumptions on the potential q(x); so we omit the proofs. 

1. If y is any solution of the SL-equation, then ( (y') 2 , —2yy, y 2 ) T is a solution of the first order system 

OH. 



(3.1) 



2. If we let a fundamental system of the Sturm-Liouville equation (|2.4|) be defined by the initial conditions 
at any xq > 

u(x 0) X) v(x 0) X) 
u'(x ,X) v'(x ,\) 

then a corresponding fundamental system of solutions of equation (|1.6p is 





" 1 


" 
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U=[U 1 ,U 2 ,U 3 ] = 



{u') 2 u'v' (v') 2 

—2uu' — \u'v + uv'] —2vv' 



(3.2) 



3. If {4>(x, A), 9(x, A)} are the Frobcnius solutions defined in Section 2 (in all the cases) and 
normalized by (I2.16[) so as to ensure that W x ((/)(■, A), 9(-, A)) = 1, then a corresponding fundamental 
system of solutions of equation (jl.6|) is 



U=[U 1) U 2 ,U 3 ] = 



{9') 2 6 

-299' -[(9'cj 
a2 



4. An indefinite inner product on the solution space of equation (|1.6p may be defined by 

(Ui,U 2 ) := 2(PiR 2 + P 2 R\) - QiQ 2 = const, independent of x e [0, oo) 
where U k = (P k , Q k , R k ), k = 1, 2. 

5. For any solution U = (P,Q,R) T of equation (|1.6[) , 

i.e. 

APR — Q 2 = const, independent of x e [0, oo) 

6. If Ui and U 2 are any two solutions of equation ()1.6[) represented in the form, 



(3.3) 



(3.4) 





' Pj ' 




- {6'f ' 




e'(j>' 




\ m 2 I 




Qj 


= a 3 


-299' 


+ bj 




+ Cj 


-2^' 




Rj 




9 2 




9<f) 







in terms of the Frobenius solutions {</>(•, A), 9(-, A)} of section 2, then we have 

(U u U 2 ) := 2(P 1 R 2 + P 2 R 1 ) - QiQ 2 - 2{a x c 2 + c x a 2 ) - b x b 2 

In particular, 

{U u Ui) := APiRi -Qf = 4a lCl - b\. 
7. Similarly, if these same solutions, U\ and U 2 , of equation ()1.6j) are represented in the form, 







(u') 2 






= d i 


-2uu' 




. R j - 




u 2 





[uv + u v\ 
III' 





r w) 2 i 


+ c j 


-2vv' 








v 2 



(3.5) 
(3.6) 

(3.7) 
(3.8) 

(3.9) 
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in terms of the solutions denned in (|3 - If) - f|3 . 2|) we have 

(U U U 2 ) := 2{P 1 R 2 + P 2 R 1 ) - Q X Q 2 = 2{a x c 2 + cio a ) - hh- (3.10) 

In particular, 

(U u U x ) := ±P X R X - Qf = 4a!C! - b\. (3.11) 
It follows from (|3~7)l . (|3TT0|) and ((3~8| . (pTTTj) that we must also have 

2(aiC2 + C1O2) — 6162 = 2(aiC2 + £1(22) — &162, and4aici — b\ = AaiC\ — b\. (3.12) 

8. If y is any solution of the SL equation (jl.ll) and t/=(P,Q,R) T is any solution of the companion 
system (|1.6j) then 

■fj-[Py 2 + Qyy' + R(y') 2 } = 0, (3.13) 

i.e., 

\)y 2 (x 1 A) + Q(a;, X)y(x, \)y'(x, A) + X)(y'(x, A) 2 = constant, independent of x. (3.14) 

4 A Spectral Density Function Characterization of Al-Naggar 
and Pearson 

In this section we consider the Sturm-Liouville problem 

- y"( x ) + f ® + - + £ J j/(x) = Ay(x), x e [A, 00), A > 0, (4.1) 

y(A) = 0. (4.2) 

We shall assume that Assumption 1, Case I, holds, so that the above potential q is continuous in 
(0,oo) and also has a continuous derivative. In addition we make the following assumption: 

Assumption 3: Near x = 00 : 

For So > we have either 

?£Li(i ,oo) (4.3) 

or 

q' € Li(xq, 00), q € ACi oc [xq, 00), and lim q(x) = 0. (4-4) 



Under the assumption (|4.3j) or the assumption (|4.4|) it was established in Fulton, Pearson and Pruess 
[T2l Thml and Cor 2] that the initial value problem (|1.6[) - (|1.7p has a unique solution for all A € (0, 00). 
Henceforth we denote this unique solution by 



Ux(x,\) 



Pi Or, A) 

Qi(x,\) 
Ri{x,X) 



(4.5) 



that is, U% is the unique solution (under Assumption 3) for which 



Pi (00, A) ^ 






Qi(oo, A) 


-(i) 


for A € (0, 00) 


i?i(oo,A) y 







lim J7i(a;, A) = 

a:— ^cxd 

The assumptions (j4.3[) and (|4.4p were used in H3] to obtain spectral density function characteri- 
zations of the type (|1.5[) when the left endpoint is regular. Al-Naggar and Pearson [TJ [2] also obtained a 
spectral density function characterization of this type when the left endpoint is regular, using a different 
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approach. Their approach was based on the determination of intervals of a.c. spectrum by locating 
those intervals of the real A- axis where subordinate solutions do not exist. The purpose of this section is 
to present their analysis as it applies to the problem (|4.1[) - (|4.2j) and under the additional Assumption 
3 on the half line [A, oo), A > 0, and show that it also guarantees a.c. spectrum for A € (0, oo) and 
yields a spectral density function formula of the type (|1.5p (see (|4.30|) below and Cor 4]); this route 
to the spectral density function characterization represents an alternative to the analysis in [12l [13] ■ As 
we shall see, the approach of Al-Naggar and Pearson has a major advantage in that it extends nicely 
to obtain a corresponding spectral density function characterization of the type (|1.5[) for the doubly 
singular equation (|4.ip on (0, oo). This will be done in the next section. 

In [TJ [2] the analysis is focused on the third order ordinary differential equation (see (14.131) below) 
satisfied by the third component of a solution of Appell's first order system; here, we modify the approach 
slightly so as to focus attention on the system (|1.6|) . so that we can properly exploit the results from 
[12] on uniqueness of the above solution Ui satisfying the initial condition (|1.7p . 

Letting {u(-, A), «(•, A)} be the fundamental system of solutions of (|4.1[) defined by the initial condi- 
tions (|3.ip at xq = A, the Titchmarsh-Weyl m-function associated with the problem (|4.ip - (|4.2p is defined 
by 

V A (-, A) := u(; A) + m A (A)u(-, A) <= L 2 (A,oc), for all 7m (A) + 0. (4.6) 
Then, as is well known, this m-function is a Nevanlinna function and therefore admits the representation 



mA{z) = a + f3z + 



1 



1 + f 



dpA(t), ae (—00,00), f3 > 0, 



where the inversion integral for pa in terms of m is the Titchmarsh-Kodaira formula 



(4.7) 



Pa(X) =\im — / Im\mA(t + ie)]dt. 



(4.8) 



Since the Assumption 3 ensures a.c. spectrum for A G (0, 00) we may differentiate (|4.8p to obtain the 
spectral density function as 



/a (A) := PaW = nm —Im[m A (\ + ie)] 

ej.0 7T 



(4.9) 



We now proceed to transform (|4.9p into the form (|1.5p ; after some analysis this yields (14.301) in 
Theorem 1 below. In the sequel it will be helpful to make use of the fundamental system of solutions of 
Appell's system (|1.6p given in Q3.2p in terms of the solutions {u(-,X),v(-,X)} of equation (|4.1|) fixed by 
the initial conditions (|3.ip with xq = A 



Lemma 1. Assume the potential q satisfies Assumption 3. 

(i) For all A € (0, 00) let U\ be be the unique solution defined at x = 00 in (|4.5[) . Then using the 
indefinite inner product on the solution space of (|1.6p defined in (I3.7P we have 

(E7i, Ui) = 4Pi(x, X)Ri(x, A) - (Q^x, A)) 2 = 4. (4.10) 

(ii) If the solution U\ is represented in the form (|3.9p . say 



(4.11) 





' Pi(x,X) ' 




- (u'f ~ 




uV 




- (v'f ' 


Ux{x,X) = 


Qx(x,X) 


= a 


-2uu' 


+ b 


— [uv' + u'v] 


+ c 


-2vv' 




_ Ri(x,X) _ 




u 2 




uv 




v 2 



then 

Adb- (cf = 4. 

PROOF, (i) To get the constant 4 observe from (|3.5p and the initial condition (|1.7I) that 

4Pi(x,A)JJi(a!,A)- (Qi(x,X)) 2 = lim 4Pi (x, X)Ri (x, A) — (Q 1 (x 7 X)) 2 = 4VX ■ (4=) 

x^oc VA 



(4.12) 



= 4. 
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(ii) The conversion of the indefinite inner product in terms of the coefficients {a, b, c} is the property 
(|3.7j) and ()3.8j) : this is readily proved by substitution of the components of U\ into the inner product 
formula and use of the wronskian relation W x (u(-, A), v(-, A)) = 1. I 

It can be shown that the third component, R(x, A), of a solution of Appell's system (II.6[) satisfies 
the third order equation of [2J p. 6584, Equa. 5], 



d R ... , ,.dR , 
_ +4{ X- q (x))—-2q'(x)R = 0. 



(4.13) 



Remark In Appell's paper this third order equation is [3j p. 213, Equa (5)]. 

Also, a general solution of (|4.13l) can be obtained from a suitably normalized solution Ri(x, A) as in 
p. 6587, Lemma 2]. Here we generalize this technique to obtain from the given solution U\, two other 
linearly independent solutions. This is the content of the following lemma. 

Lemma 2. Let U\ = (Pi,Qi,Ri) T be the unique solution of ()1.6[) defined at x — oo by the initial 
condition (|1.7|l . Then we can write the general solution of (|1.6p in the form 



U = 



P 

Q 

R 



= faU 1 +faU 2 +faU 3 



where 



U 2 (x,X) 



U 3 (x,\) 



and 



PiCOs2 7 + (Qi/i?i)sin27- (2/Pi)cos27 
Qi cos 27 + 2 sin 27 
R\ cos 2 7 

Pi sin 27- (Qi/i?i)cos27 - (2/Pi)sin27 
Qi sin 27 — 2 cos 27 
i?i sin 27 



j(x) = / l/R 1 (t,X)dt 

Jx 



(4.14) 

(4.15) 
(4.16) 



for some xq > 0. 
Proof. The fact that 

R(x, A) := Ri {x, A) + fa cos 2 7 + fa sin 2 7 ] 

is a general solution of (|4. 13[) follows as in [2, Lemma 2, p. 6587]. Since (|4.13l) is satisfied by the 
third component of any solution of (|1.6p . we can generate a general solution for (|1.6[) by computing 
Q(x, A) = —dR/dx and then P(x, A) = [— dQ/dx + 2(A — q)R}/2 and expressing these solutions of (I1.6P 
in terms of Pi, Qi, R\. This gives the result (|4.14p . Alternatively, a direct substitution of U 2 and U3 
into (|1.6p and use of the formulas for P[, Q[ and R[ will verify that U 2 and U3 are solutions of (|1.6p . 



Corollary, (i) The general solution (|4.14p in Lemma 2 satisfies 

(U, U) = APR -Q 2 = 4((/3i) 2 - (fa) 2 - (fa) 2 ). 

Similarly, if 

U = (P, Q, R) = faUi + faU 2 + faU 3 , 
we have for the inner product that 

(U, U) = 2(PR + RP) - QQ 

= A(fafa-fafa-fafa). 



(4.17) 



(4.18) 



(ii) The solutions {U\, U 2 , U3} in Lemma 1 are mutually orthogonal with respect to the indefinite inner 
product defined in (I3.4p . 
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PROOF, (i) For (|4.17j) calculate APR-Q 2 using (|4.14j) and the normalization (I4.10p . This simplifies to 

4(/3i +^ 2 cos27 + /3 3 sin27) 2 - 8(^2 cos2 7 + j3 3 sin2 7 )(/3 1 + 2 cos2 7 + #5 sin27) 

- 4((/3 2 ) 2 sin 2 2 7 + (/3 3 ) 2 cos 2 2 7 - 2/3 2 /3 3 sin 2 7 cos 2 7 ) 
= 4((/3 1 ) 2 -(/3 2 ) 2 -(/?3) 2 ). 

the proof of (|4.18p is similar. 

(ii) Taking U = U\, which statisfies the normalization (|4.10l) . and U = U 2 in (|4.15p we have from (|4. 1 8[) 
that (Ui, U 2 ) = 4(1 • - • 1 - • 0) = 0. Similarly, (Ui,U 3 ) and (U 2 , U 3 ) are zero. ■ 

Remark. The third component, Ri, of (|4.1ip satisfies the third order equation (|4.13|) and since 
R\{x, A) — > 1/a/A > 0, it is in fact the same quadratic form as Y(x,X) which was employed in [2 
Lemma3 and Thm2]. 

The following lemma gives a limit relation involving the solution, U± , and the solutions from Lemma 
2 which are orthogonal to it. 

Lemma 3. We assume the potential q satisfies Assumption 3. 

(i) For the solution U\ of the initial value problem (|1.6[) - (|1.7p . let U 2 and U 3 be the linearly independent 
solutions of Lemma 1 which are generated by using (Pi, Q\, R\) T in (|4.15p and (|4.16p . Then, for any 
linear combination 

U= ( Q j = (3 2 U 2 + f3 3 U 3 

we have for all xq > 0, 

f N R(x. X) dx 

U m =0. (4.19) 

(ii) The representation (14. lip for U\ has 5(A) > and c(A) > for all A G (0, 00). 

PROOF (i). The solutions U 2 and U 3 in (|4.15[) and (|4.16l) have R 2 = Ri(x, A) cos(2 7 (a;)) and 
^3 = Ri{x, A) sin(2 7 (x)) , where 

j(x)= I [1/R l {t,\)]dt ) 

Jxq 

for all x > 0. The fact that Ri(x) > for all x > follows from (|4.20p below. Consequently, it suffices 
to prove that (|4.19p holds for these two choices of R. For R 2 we have 

r N r N d 

/ Ri cos 2 7 dx = / 0.5(i?i) 2 — sin27dx 

Jx Jx 

= 0.5( J Ri) 2 sin2 7 |^ + / R 1 Q l sin2 7 dx. 



Since Qi(a;) — > as a; — > 00 we can prove that 

f R1Q1 sin 27 da; 

Km Jao ^ — = 0. 

N ^ f N Ri dx 

Jxq 1 

Given e > 0, pick x t sufficiently large that |Qi(x)| < e for x > x e . Then we have 



N 

R1Q1 sin 27 dx 

x 



■ N 

< M I Ridx + e I Rx dx, 



where M is a bound on Qi in [^o^e]- Now pick N sufficiently large that 



M Ridx 

J Xg 

N 

Ri dx 

Xo 



< e, 
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which is possible since f Ri dx = oo. Then the above quotient is less than 2e. For the boundary term 
we have 



[i?i(iV,A)] 2 = [i?i(iro,A)] 2 + / 2R 1 R' 1 dx 

J x a 

= [i?i(a;o,A)] 2 - / 2R 1 Q 1 dx, 

J In 



SO 

Urn lR f A)? =»- 
N ^°° [ N R x dx 

Jx 1 

by employing the same argument. The proof for R3 is similar, 
(ii) In the representation (|4. 1 1|) for JJ\ we have 

Ri(x, A) = au(x, A) 2 + bu(x, X)v(x, A) + cv(x, A) 2 . 

To see that a > and c > for all A <E (0, 00) we first observe from (|4.12[) that 45c — (b) 2 = 4 requires 5 
and c to be of the same sign; and this must hold for all A £ (0, 00) since they are continuous and cannot 
pass through zero (which would violate (|4.12l) '). Factoring out 5 we have 

Ri(x, A) = a[u(x, A) 2 + (b/d)u(x, X)v(x, A) + (c/a)v(x, A) 2 ] 

with the coefficients of u 2 and v 2 positive. Therefore R\ admits a factorization of the form 

Ri(x, A) = d[(u(x, A) + av(x, X))(u(x, A) + av(x, A))] 

= a\u(x, A) + av(x, A)| 2 , (4.20) 

where a = a\ + iui must satisfy (by (|4.12[) ) 

01 = -6/(25) and a\ = (c/~a) - [(6) 2 /(4(5) 2 )] = l/(a) 2 . (4.21) 

Since a, b, and c are real, we must have either U2 = 1/5 or ai = —1/5; but either way the factorization 
remains the same with a and a switched. Since \u — av\ 2 > 0, it follows from the above factorization 
that 5(A) > for all A £ (0, 00); otherwise, lim :c _ i . 00 Ri(x, A) < contradicting the fact that i?i(oo, A) = 
l/VA > 0. Hence, also 5(A) > for all A £ (0, 00). I 

Since the definition of 012 must be ±1/5(A) from (|4.21l) and since this indeterminacy is actually 
immaterial, we choose to take 012 = 1/5(A), so that 

a(X):=~+il (4.22) 
la a 

Using the fundamental system {u(-, A), v(-, A)} we now define the complex-valued solution of (|4.ip for 
real A £ (0, 00), 

if>A(x, A) := u(x, A) + a(X)v(x, A). (4.23) 

The key idea of Al-Naggar and Pearson, which enables identification of a.c. spectrum, is embodied in 
the following requirement: 

Definition. The general Sturm-Liouville equation (11.11) satisfies Condition A, for a given real value of 
A, if and only if there exists a complex-valued solution y{x, A) of (|1.1[) for which 

N 

2 



y(x, A) dx 

lim = for x > 0. (4.24) 



N -y 00 



\y(x,X)\ dx 

x 
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We now prove that equation (|4.1[) satisfies Condition A for all A G (0, oo). 



Lemma 4. Assume q(x) in (|4.1j) satisfies Assumptions 3. Then for xo > and all A £ (0, oo) 



r/u(a;, A) 2 dx 

l im ; =0 . (4.25) 



£ \Mx,\)\ 2 dx 



Proof: Since equation (|4.13l) is satished by all linear combinations of v(-, A) 2 , u(-, X)v(-, A) and u(-. A) 2 
(see [H Lemma 1, p. 6584]), it follows that Vm(',A) 2 is a solution of (I4.13|) and since it is complex- 
valued, also that Re [4>a(-, A) 2 ] and Im [iPa(', A) 2 ] satisfy (|4. 13|) . Accordingly, it follows that there exist 
solutions U2 = (P2, Q2, ^{^Pa)) T an d U3 — (P3, Q3, lm(tp\)) T of the first order system (|1.6[> . since any 
real solution i? of (|4.13l) can be used to generate a solution of (|1.6I) having i? as its third component; 
e.g., let Q' = -R! and P = -Q'/2 - (A - q)R. From and (|4~2U1) we readily deduce the following 

representations of i?i(-, A), i2e[^(-, A)] 2 , and Jm[Vu( - , A)] 2 of the form (|4.11l) (or, the third component 
of KHl ): 

Ri(x, A) = au 2 + 2aa\uv + a(a 2 + a 2 )w 2 
Re [ipA(x, A)] 2 = it 2 + 2ai?iw + (a 2 - a|)u 2 
Imfy^x, A)] 2 = 2ct2uv + 2ai(X2V 2 . 



It now follows from these formulas that the above solutions U2 and C/3 associated with Re[ipA(-, A)] 2 and 
Im[ipA('i A))] 2 are orthogonal to U± in the sense of the inner product defined in Q3.7[) . i.e., we have 

(Ui,U 2 ) = 2[5(a 2 - o|) + a(a 2 + o|)] - 4aa 2 = 0, 

and 

(Ui,U 3 ) = 2[5(2aia 2 ) + 0] - 45aia 2 = 0. 
Hence it follows from Lemma 3(i) that for all x a > 

/^Re[^(x,A)] 2 ^ 
lim — i — jr? = 

and 

,. J lm *Pa(x, X)] 2 dx 
hm T7 = 0, 

from which (|4.25[) follows. I 

Assumming the Sturm-Liouville equation has a potential g which is LP at x = 00 and regular 
at a finite left endpoint, Al-Naggar and Pearson obtain the following results in p] Lemmal and Thm2] 
(where the fundamental system {u, v} is defined by initial conditions at the left endpoint so that v 
satisfies a general regular boundary condition and W x (u(- 1 A), v(-, A) = 1): 

Lemma 5 (Al-Naggar and Pearson). Let Id (—00, 00) be an interval on which Condition A holds for the 
general equation 11.1% and let the fundamental system {u(-, A), v(-, A)} be defined by the initial conditions 
nn\) at x = A. Then 

(i) There exists a complex valued function M(X) on I which is uniquely defined by the properties: 

N 

2 



(u(x,X) + M{X)v(x,X)Y dx. 
(a) Im[M(X)] > and (6) Jim ^ = 0. 



AT-s-oo 



\{u(x,X) + M(X)v{x,X)\ 2 dx 

x 
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(ii) For A £ I the function M in (i) is the boundary value of the Titchmarsh- Weyl m-fuction defined by 
&4-6\ ), that is, 

M(X) = lim[77u(A + ie)]. 

Proof. Statements (i) and (ii) are, respectively, Lemma 1 and Theorem 2 from pQ. B 

We now apply these results to the problem (I4.1l) - (|4.2p . 

Lemma 6. Assume that for the problem (|4.ip - (|4.2l) Assumption 3 holds. Then with a(A) defined by 
(|4~2"2"|) we have for all A G (0, oo) 

a(X) =1im[m A (\ + ie)]. ■ (4.26) 

Proof. By the uniqueness of M(A), and the fact for all A G / = (0, oo) a(X) has positive imaginary 
part (see (4.22) and Lemma 3(h)) and ij} A (x, A) satisfies property (b) in Lemma 5(i) (see Lemma 4) it 
follows from Lemma 5(h) with / = (0, oo) that for all A g (0,oo) we must have (|4.26[) . I 



We are now ready to prove the representation of f a (A) in the form (|1.5|) . 

Theorem 1 Assume the potential q is given as in (|4.1j) and that Assumption 3 holds. Let a(A) be 
dchncd as in (|4.22p and iPa(', A) as in (|4.23l) . Then the spectral function defined by (|4.8[) for the problem 
(|4.1[) - (|4.2[) is absolutely continuous for A G (0, oo) and the corresponding spectral density function admits 
the following representations for A G (0,oo): 

/ A (A) := p'(X) = lim -Im[m(X + ie)} (4.27) 

= «?W (4.28) 

' (4.29) 



7ra(A) 

~ 7r[Pi(ar, X)(v(x, A)) 2 + QJx, X)v(x, \)v'(x, A) + i?i(x, X)(v'(x, A)) 2 ] 

Proof: The statement (|4.28p . and the absolute continuity of the spectral function p, follows from 
Lemma 6 by taking imaginary parts on each side of equation (|4.26p . The statement (I4.29P follows from 
the definition of ct2 in (|4.22|1 . To obtain (|4.30p substitute U\ = (Pi,Q±, R\) T from the representation 
(14. lip in terms of {u(x, A), v(x, A)} and collect coefficients of a(A), 6(A) and c(A) to obtain 

Pv 2 + Qvv' + R{v') 2 = a{X)[W x (v,u)] 2 = a(X). 



Remark. Putting x = in (|4~30| yields f A (X) = 1/(tt.Ri(A, A)) which was a main result of Theorem 
2 in [2]. 

The following well known spectral density function formula, due to Titchmarsh [23], 1946, and Weyl 
[24], 1910, in the case of the assumption (|4.3p (and due to Pearson [19l [H] in the case of assumption 
(|4.4p ) also follows readily from Theorem 1: 

Corollary. Under the assumptions of Theorem 1 we have for all A G (0, oo) 

f A (X) = lim = — ■ . (4.31) 

x^n[VX(v(x,X))2 + j=(v>(x,X))2] 



Proof: The Assumption 3 ensures (see [12] Thm 2]) that the solutions v and u defined by the initial 
conditions (|3.1I) are bounded for sufficiently large x. Hence, making use of the initial condition (I1.7P 
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which U± satisfies, we have 

Pi(x, X)(v(x, A)) 2 + Qi(x, X)v(x, X)v'(x, A) + Ri(x, X)(v'(x, A)) 2 
= lira [Pi(x, X)(v(x, A)) 2 +Q 1 {x,X)v(x,X)v'(x,X) + R 1 {x,X){v'{x,X)) 2 } 

x— >oo 

= lim Q 1 (x,X)v(x,X)v'(x,X)+ lim [Pi (a:, X)(v(x, A)) 2 + Ri(x, X)(v'(x, A)) 2 ] 
= + lim [V\(v(x, A)) 2 + -^(v'(x, A )) 2 ]^ 

so it follows that (|4.30[) gives rise to the characterization ()4.31|) . I 

Remark. In [12] we made use of the formula (|4.31j) to establish (|4.30|) . Here, by linking the spectral 
density function first to the m-function, and following the approach of Al-Naggar and Pearson, we have 
obtained a direct proof of (|4.30j) . from which the older result (14.311) follows as a consequence. 

5 Generalization of a Spectral Density Function Characteriza- 
tion to Doubly Singular problems 

In this section we consider the Sturm-Liouville problem 

ry := -y"(x) + 1 1| + | + £ q n+2 x n J y(x) = Xy(x), x € (0, oo), (5.1) 
W ((y(;X),<P(;0)) = limW x ((y(;X),(P(;0)) = 0, ifx = 0isLC, (5.2) 

x— ^0 

where </>(x, 0) is the first Frobenius solution for A = given in (|2.16j) . Since this is also the principal 
solution at x = in all (LC and LP) cases, the boundary condition (I5.2p is the Friedrichs boundary 
condition in the LC cases at x = and selects 4>(x, X) for all A G C; this boundary condition is also 
automatically satisfied by (f>(x, A) in all the LP cases at x = 0. In this section we adopt the Assumption 
1 and 2 from Section 2, and the Assumption 3 from Section 4. Our aim is to extend the spectral 
density function characterization in Theorem 1 (under the above 3 assumptions) to obtain the formula 
(|1.5[) for the doubly singular problem (I5.1[) - (|5.2p . The Assumption 1 ensures that x = is a singular 
point of type LP/N or LC/N; the Assumption 2 ensures that x = oo is a singular point of type 
LP/O-N with cutoff value A = 0; and the Assumption 3 ensures that we have a.c. spectrum in 
(0,oo). The underlying self-adjoint operator associated with (|5.ip - (|5.2|) has the domain 

D(A);=(feL 2 (0,oo) \ f(x) e AC/ oc (0, oo), r/eL 2 (0,oo), 

Hm^ o W^(/(-) 5 0(->O))=o| (5.3) 

in the LC cases at x = 0, and 

D(A):=|/€L 2 (0,oo) | / e AC? oc (0, oo), r/eL 2 (0,oo)| 



(5.4) 



in the LP cases at x — 0. The associated eigenfunction expansion theory which obtains the eigenfunction 
expansion in the form (|1.2|l for both of the above cases was given in [10] and [14] , and explicit formulas 
for the corresponding Titchmarsh-Weyl m-function and associated scalar spectral function were obtained 
in [Tl] for all cases of the special potential 

q(x) = % + — , qo and q x satisfying (I2.2I) . (5.5) 
x 2 - x 

Here the Titchmarsh-Weyl m-function is defined as in [101 [Tl] by 

*(-,A) := 0(-,A) -m(A)0(-,A) € L 2 (x ,oo), x a > 0, for all Jm(A) ^ 0. (5.6) 
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where (/>(■, A) and 6(-, A) are the first and second Frobenius solutions normalized as in (|2.16[) . We now 
repeat some basic information from (141 concerning the Titchmarsh-Weyl m- functions defined by (|5.6|) . 
In the LC cases at x — the m-function defined by ()5.6j) is a Nevanlinna function (that is, a function of 
class No), and therefore admits an integral representation of the form (|4.7|) . and the inversion integral 
for p in terms of m is a Titchmarsh-Kodaira formula like (|4.8[) . The eigenfunction expansion for the 
problem (|5.ip - (|5.2l) when x = is LC has the form (| 1 . 2|) where </>(•, A) is the first Frobenius solution in 
(|2.16p . and p is the spectral function obtained from the above m-function. 

In the LP cases at x = 0, the m-function defined by (|5.6I) is a generalized Nevanlinna function of 
class N K for some k > 1 (see [HI p. 188]) and it follows using the theory of Krein and Langer [TH] for 
these functions (see [H Thm 3.5 and Lemma 4.1]) that they admit the representation 



o 



where n,m> 1, ay £ {— oo, oo), a m ^ if m > 0, and where a is a measure on R satisfying 

da(t) 



1 + i 2 



< oo. 



The spectral function for the associated self-adjoint operator A in (|5.4[) is then defined interms of a as 

r-A 

2\ n 



p(X):= (l + s 2 )da(s), AG (-oo.oo). (5.7) 

While the Titchmarsh-Kodaira formula is well known for No functions, it was only recently established 
in various cases with two LP endpoints and simple spectrum by Gesztesy and Zinchencko [16] and 
Fulton and Langer [TJj. For the case of equation (|5.1I) . where generalized Nevanlinna functions of class 
N K arise, we quote this result, and also relate it to the classical real-variable definition of Levitan and 
Levinson (see [H Thm 4.7 and 4.8]): 

Theorem 2 (Fulton and Langer). Consider equation (15. ip with x — of LP type and suppose Assump- 
tions 1 and 2 hold. If A, Ao are points of p-measure zero (not discrete eigenvalues of the associated 
self-adjoint operator A), then the spectral function defined by (|5.7p has the representation 



1 f x 

p(X) — p(Ao) = lim — / Im\m(p + it 

^° 71" J\ Q 



)}dp (5.8) 



lim y 



\ ]b e{\ ,\)r\<j(A b ) Jo \(t>( x ^jb)\ 2 dx 



Note: When (|5.ip is of LC type at x — we have n = in (|5.7p and a standard Nevanlinna representation 
of m(A), for which the inversion formula is also (|5.8p . 



Here, At is the corresponding truncated self-adjoint operator on (0, b] with any regular boundary 
condition at x = b. It follows from [Mj Thm 4.5] that for the LP cases at x = the function defined 
by (|5.7p or (|5.8p is the spectral function which arises in the eigenfunction expansion (|1.2p . For proofs of 
convergence results and Parseval relation we refer to (TOl H3] . 



Remark. The above theorem justifies the use of the second formula in (|5.8p for the computation of the 
spectral function p as was implemented in the software package SLEDGE when both endpoints are of 
LP type. For Sturm-Liouville problems satisfying Assumptions 1,2 of Section 2, SLEDGE [15, Equa 
(1.13)] does in fact normalize the (^-solution as in (|2.6p . (|2.9p and (|2.12l) . 

If we make the Assumption 3 in addition to the assumptions 1 and 2 of Section 2, the spectrum 
is a.c. on (0,oo). In the case when the left endpoint is regular, it was shown under Assumption 3 in 
P2J Thm 1 and Cor 2] that the initial value problem (|1.6p ~ (|1.7p at x = oo has a unique solution for each 
of the cases (|5.3p and l|5.4p : also, that the spectral density function has the form (|1.5p (see [HI Cor 4] 
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and [13l Thm 1]) for each of these cases. The proof required linking the formula (|1.5p to the well known 
result of Weyl and Titchmarsh (|4.31[) in the Corollary to Theorem 1. An alternative approach to the 
proof of (|1.5|) when the left endpoint is regular was given by Al-Naggar and Pearson p] [2] as described 
in Section 4. We now follow this method of analysis to generalize the above Theorem 1 to the case when 
both endpoints are singular and all three assumptions hold. Since the spectral function is absolutely 
continuous in (0, oo), we may differentiate in (|5.8I) to obtain the spectral density function as 

/(A) := p'(X) = lim -Im[m(X + it)]. (5.9) 

The eigenfunction expansion associated with the underlying self-adjoint operator A when x = is LP 
has the form (|1.2I) where (/>(■, A) is the first Frobenius solution in (|2.16[) . and p is the spectral function 
defined in (|5.7|l or (j5.8[) . In both cases (LC and LP at x = 0) the spectral density function is given 
by (|5.9[) in terms of the m-function which is defined in (|5.6I) above relative to the suitably normalized 
Frobenius solutions. 

We proceed now to transform (|5.9p to the form (|1.5[) ; after some analysis this yields (|5.30l) in Theorem 
3 below. Since many of the necessary lemmas which are required are the same as in Section 4, we list 
and prove only those lemmas which undergo some change as a result of allowing the left endpoint x — 
to be a singular endpoint satisfying Assumptions 1,2. We also adopt the notational convention that 
Lemma n* in this section represents the analog of Lemma n in Section 4. 

In the sequel it will be helpful to make use of the fundamental system of solutions of Appcll's system 
(|1.6p given in (|3.3[) in terms of the suitably normalized Frobenius solutions {${■, A), &(■, A)} of equation 
(52) fixed by the definition ||23B)> . 

Lemma 1*. Assume that for equation (|5.1[) Assumption 1, Casel, and Assumption 2,3 hold. 

(i) Same as Lemma l(i); this yields 

(E7i,[/i) =4P 1 (x 1 \)R 1 {x,\)-(Qi{x 1 \)) 2 =4. (5.10) 

(ii) If the solution U\ is represented in the form (|3.6|) . say 





' Pi{x,X) ' 




- (er ' 




&4>' 




r m 2 i 


Ui(x,X) = 


Qi(x, A) 


* 

= a 


-299' 


+ b* 


-[eft + e'<f>] 


+ c* 


-2j><f>' 




. Ri(x,X) _ 




9 2 




9(j) 







(5.11) 



then 

4a*c* - (6*) 2 = 4. (5.12) 

Proof, (ii) The conversion of the indefinite inner product in terms of the coefficients {a*, b* , c*} is the 
property ()3.7p and (|3.8j) : this is readily proved by substitution of the components of U\ into the inner 
product formula and use of the wronskian relation W x (<p(-, A), #(-, A)) = 1. I 

Lemma 3*. Assume that for equation (|5.1j) Assumption 1, Casel, and Assumption 2,3 hold. 

(i) Same as Lemma 3(i). 

(ii) The representation (15. 1 1[) for U\ has a* (A) > and c*(A) > for all A e (0, oo). 
Proof: (ii) In the representation (|5.1ip for U\ we have 

Ri(x, A) = a*9{x, A) 2 + b*6(x, \)<j>(x, A) + c*4>{x, A) 2 . 

To see that a* > and c* > for all A e (0, oo) we first observe from (|5. 12[) that 4a* c* - (6*) 2 = 4 
requires a* and c* to be of the same sign; and this must hold for all A £ (0, oo) since they are continuous 
and cannot pass through zero (which would violate (|5.12l) ). Factoring out a* we have 

Ri(x,X) = a*[9(x,X) 2 + (b*/a*)6(x,X)(p(x,X) + (c*/a*)<fi(x, A) 2 ] 

with the coefficients of 9 2 and </> 2 positive. Therefore Ri admits a factorization of the form 

R^x, A) = a*[(9(x, A) - £<j>(x, X))(6(x, A) - &(x, A))] 

= a*\6(x,X) -£<f>(x,X)\ 2 , (5.13) 
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where £ = £1 + i& must satisfy (by (|5.12|l ) 



& = -6*/(2a*) and ^ - (c*/a*) - [(6*)V(4(o*) 2 )] = VK)- (5.14) 

Since a*, 6*, and c* are real, we must have either £2 = l/o* or £2 = — l/ a *i but either way the 
factorization remains the same with £ and £ switched. Since \9 — £</>| 2 > 0, it follows from the above 
factorization that a* (A) > for all A € (0,oo); otherwise, lim :E _j. 00 Ri(x, A) < contradicting the fact 
that i?i(oo, A) = 1/vA > 0. Hence, also c*(A) > for all A £ (0,oo). I 

Since the definition of £2 must be ±l/a*(A) from (|5.14j) and since this indeterminacy is actually 
immaterial, we choose to take £2 = 1/a* (A) , so that 

£(A) + . (5.15) 

sv ; 2a* a* v ; 

Using the Frobenius solutions normalized by (j2.16[) we now define the complex- valued solution of (I5.1[) 
for real A £ (0, 00), by 

ip(x, A) := 6{x, A) - i{\)(f){x, A). (5.16) 

We now prove the following lemma that equation (|5.1[) satisfies Condition A for all A £ (0, 00). 

Lemma 4*. Assume q(x) satisfies Assumption 1, Case I, and Assumptions 2, 3. Then for Xo > 
and all A £ (0, 00) 

f i/)(a;, A) 2 

Um J *° ^ ; =0. (5.17) 

^°°£|^,A)|2& 

PROOF: Since equation (|4.13p is satisfied by all linear combinations of </>(•, A) 2 , #(•, A)</>(-, A) and #(•, A) 2 
(see [H Lemma 1, p. 6584]), it follows that ip(-, A) 2 is a solution of (|4.13[) and since it is complex- 
valued, also that Re [^(-,A) 2 ] and Im A) 2 ] satisfy (|4.13p . Accordingly, it follows that there exist 
solutions U2 = {P2, Q2, Re(ip 2 )) T and U3 = (P3, Q3, Im(i/' 2 )) T of the first order system (|1.6|) . since any 
real solution R of (14.13|) can be used to generate a solution of (|1.6|) having R as its third component; 
e.g., let Q' = —R' and P = —Q'/2 — (A — q)R. From (|5.13[) and (|5 . 1 6[) we readily deduce the following 
representations of i?i(-, A), Re[ip(-, A)] 2 , and Im[tp(-, A)] 2 of the form (|5.11[) (or, the third component of 

JEHU): 

,* Q 2 o„*<- dX 1 „*lc2 1 £2w2 



Ri(x, A) = a*8 z - 2a*ii6<l> + o*(£f + 

2 2 ) 



Re [</>(x, A)] 2 = # 2 - 2&00 + Hi - £* v " 



In#(x,A)] 2 = -2^ + 2& 6<^ 2 



It now follows from these formulas that the above solutions U2 and U3 associated with Re[ip(-, A)] 2 and 
Im[ip(-, A))] 2 are orthogonal to U\ in the sense of the inner product defined in (|3. 7[) . i.e., we have 

{U U U 2 ) = 2[a*(e i - CI) + a*(Zi + $)] - 4a*g = 0, 

and 

(U U U 3 ) = 2[a*(2&6) + 0] - 4a* = 0. 
Hence it follows from Lemma 3(i) that for all Xo > 

J x N o Remx,X)] 2 dx 
km — ^ = 



and 



lim 



J xo Im[iP(x,\)} 2 dx 



f i?i(x,A)da; 



17 



from which (|5 . 1 T[) follows. H 

Lemma 5* (Al-Naggar and Pearson). We assume for equation &5.1\) that Assumption 1, Case I, and 
Assumptions 2,3 hold. Let I C (— oo, oo) be an interval on which Condition A holds. Then 
(i) There exists a complex valued function M(X) on I which is uniquely defined by the properties: 

N 

2 



(8{x, A) — M(X)(j)(x, A)) dx. 

(a) Im[M (A)] > and (b) lim = o, for all x Q > 0. 

' \9(x,X) - M{X)(j)(x,X)\ 2 dx 



where {</>(■, A), 6(-, A)} are the Frobenius solutions defined in (|2.f 61) . 

Proof: The proof is the same as that for Lemma 1 in [T]. The proof of this lemma does not depend in 
any essential way on the choice of the fundamental system of equation (15 . L [) . 8 

Unfortunately, the corresponding statement (ii) from Lemma 5 does not carry over immediately to 
the doubly singular problem (15.ip ~ (|5.2p by borrowing information from pQ; particularly, the proof of 
Theorem 2 in [1] makes use of asymptotic behaviour of solutions which are fixed by initial conditions at 
a regular left endpoint, and therefore do not apply to the Frobenius solutions {<fr,9}. The Titchmarsh- 
Weyl m- function (|5.6p was first introduced in the papers [16j [lOj [14] : and it wasn't discovered to be a 
generalized Nevanlinna function in the LP case at x = until the 2010 paper (Tl]. A direct generalization 
of Theorem 2 of pQ for cases when the left endpoint is singular remains unknown. However, we can 
recover the analogue of part (ii) of Lemma 5 by making appeal to the uniqueness result in Lemma 5*(i) 
and gleaning information on the boundary behaviour of m(z) from known information on the boundary 
behaviour of rriA{z). This is the objective of Lemma 6*. To this end, it will be helpful to introduce 
notation for the boundary values of the "regular" and "doubly singular" Titchmarsh-Weyl functions 
rriA{X) and m(A) and for the corresponding ^-functions defined in (j4.6[) and (|5.6[) . 

Definition. Associated with the m- functions, 771^4 (A) and m(A), for the problem with regular left 
endpoint and for the doubly singular problem, respectively, we define for all x € (0, 00) and all A 6 (0, 00): 

mt (A) := \imm A (X + it)), (5.18) 

*i(a:, X) := lim * j4 (A + ie) = u(x, X) + m\(X)v(x, A), (5.19) 

cj.0 

m+(A) :=limm(A + ie), (5.20) 

cj.0 

*+(x,A) :=lim*(A + ie) = 9{x, X) - m + (X)(j>(x, A), (5.21) 

ej.0 

Lemma 6*. We assume for equation (15.11) that Assumption 1, Case I, and Assumptions 2,3 hold. 
Then 

(i) for all x € (0, 00) and all z with Im z ^ 

•''-^sk)' (5 - 22) 



(ii) For all x e (0, 00) and all A € (0, 00) 

= -cf ) >(A,X)+m+(X)0(A,X) 



* + (*,A)= 1I/A f AX \ X L^ rr . (5.23) 



(iii) For xq > and all A € (0, 00) 



[ N V+(x,X) 2 dx 

lim -^S —— = 0. (5.24) 
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(iv) For all A e (0, oo) 

Im m+(A) - Imm+(A) 

+ W |m+(A)^(AA)-^(4,A)| 2 1 J 

(v) For all A 6 (0, oo) the complex valued function £(A) denned in f|5.15[) is the boundary value of the 
Titchmarsh-Weyl m-function, that is, 

£(A) = hmm(A + ie). (5.26) 

ej,0 



Proof. For (i) make use of the fact that for Imz ^ both ^a{x,z) = u(x,X) + mA(X)v(x, A) and 
^(x,z) = 6(x,X) - m(X)(f>(x > X) are in L 2 (xo,oo), and therefore linearly dependent. Using W x (<j>,8) = 
W x ((f>, \t) — 1, the relation of linear dependence is found to be 

W(ie, Z) — 



W x {4>,u + m A {v))' 



and (15 .22[) follows on evaluation of the denominator at x = A making use of the initial conditions (|3.1I) . 
For (ii) put z = A + ie in (|5.22j) and pass e — > 0. For (iii) put (|5 . 23[) into (|5.24[) and factor the constant 
terms out of the integrals to get the equivalent statement 



f x N o *\(x,Xfdx 



lim 



But ^\{x, A) is known to satisfy Condition A by Lemma 5(h), Lemma 6, and Lemma 4 with (i.e. by 
Theorem 2 of PQ applied to the problem (|4.i p -(|4.2 p ). Hence it follows that (|5.24p also holds for all 
A E (0, oo) and all xq > 0; in other words, ^> + (x, A) also satisfies Condition A. For (iv) use 

m+(X)=W x (6(;X),*+(;X)) 

and substitute the right hand side of (|5.23|) for ^ + (x,X), evaluating W x (9,u) and W x {4>, v) at x = A 
using the initial conditions p.l[) . to obtain for all A G (0, oo) 

e(A,X)m+(X)-e'(A,X) 

m (A) = 



0(AA)m+(A)-0'(AA) 

Then (|5.25|) follows by taking imaginary parts on both sides. The right hand side of (|5.25[) is positive for 
all A € (0,oo) because the numerator is positive by Theorem 1 (equation (|4.28|) ^1 . and the denominator 
never becomes zero (by separating real and imaginary parts and observing that the cases </>(A, A) = 
and <p(A, A) ^ both yield a positive denominator). To prove (v) we observe first that for all A € (0, oo), 
m + (A) has positive imaginary part (by (15.251) ) and ^ + (x, A) satisfies (by (|5 . 24[) ) property (b) in Lemma 
5*(i). Similarly, the function £ defined in (|5.15l) has positive imaginary part (see Lemma 3 (ii) ) for all 
A £ (0, oo) and the solution ip(x,X) defined in (|5.16[) satisfies (see Lemma 4*) property (b) in Lemma 
5*(i). Hence by the uniqueness of the function M(X) satisfying the two properties of Lemma 5*(i), the 
functions m + (A) and £(A) must be identical, that is, (|5.26l) holds. I 

We are now ready to prove the representation of /(A) in the form (jl.5[) . This represents the "doubly 
singular" analogue of the spectral density function characterization of Theorem 1. 

Theorem 3 We assume that Assumption 1, Case I, and Assumptions 2,3 hold. Let £(A) be 
defined as in (|5.15|) and ifj(-, A) as in (|5.16|) . Then the spectral function defined by (|5.8[) for the problem 
(|5.1[) - (|5.2[) (in both the LC and LP cases at x = 0) is absolutely continuous for A £ (0,oo) and the 
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corresponding spectral density function admits the following representations for A € (0,oo): 



/(A) := p'(X) = lim -Im[ra(A + ie)} (5.27) 

ej.0 7T 



6(A) 



7T 
1 



7ra*(A) 



tt[P x {x, \)(<p{x, A)) + Q x {x, X)<j>{x, \)<j>'{x, A) + R x {x, \)(<j>'{x, \)Y 



(5.28) 
(5.29) 

(5.30) 



Proof: The statement (|5.28[) . and the absolute continuity of the spectral function p, follows from 
Lemma 6*(v) by taking imaginary parts of each side of (15 -26[) . The statment (j5 -29[) follows from the 
definition of £2 in (|5.15j) . To obtain (|5.30|) substitute U\ — (Pi,Qi,Ri) T from the representation (|5 . 1 1|) 
in terms of {9(x, A), (f>(x, A)} and collect coefficients of a* (A), b* (A) and c*(A) to obtain 

P(0) 2 + QU' + RW) 2 = a*(\)[W x (<f>, 8)} 2 = a* (A). ■ 



6 Some Examples with Explicitly Known Spectral Density Func- 
tions 

In this section we give the explicit formulas from 10, 14 for the Frobenius solution </>(•, A), the Titchmarch- 
Weyl m- function, and the spectral density function for some examples of the special potential (|5.5[) - We 
restrict attention to those cases which will be used as test problems for the numerical algorithms in 
Sections 8 and 9. 

Example 1: [qo = —a, a > 0; qi = £(£ + 1)] Hydrogen Atom 

-y"+ f-- + ^-t!2 )y = \y a>0, < x < 00. (6.1) 



x 

The first Frobenius solution with normalization (12.91) is 



<{>(X,\) = ^ +1 [l + J>n(AK l ] 
n=l 

= x e+1 e lxV *M(£ +l-p,2l + 2, -2ixVX) 

= \=-7-rM [}i , + i{~2ixyf\), (6.2) 

with /3 := iaj2\f\ for all A S (E. Here M is the confluent hypergeometric function of first kind and M is 
the corresponding Whittaker function of first kind. The coefficients a„(A) are polynomials in A of degree 
[n/2] which are generated from the recurrence relation 

fl " (A) = -n(n+ a 2e+l) an - liX) ~ n(n+ X 2t+l) an - 2{X) > 

and the first three are 

a 

ai 
a 2 (A) 
03(A) 



2^ + 2' 

a 2 -2(l+l)A 
2!(2£+2)(2£ + 3)' 

-a 3 + (6£ + S)aX 
3!(2£+2)(2£ + 3)(2^ + 4) : 



20 



The Titchmarsh-Weyl m-function arising from (|5.6p is 

m i (X) = k e (X) -alog(-2iVx)-a^(l-ia/(2Vx))-2-/a + iVX +pt(X), (6.3) 
where "J is the psi or digamma function, 7 is Euler's constant, 



h(X) :-- 
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and where pi(X) is a polynomial of degree I (see [TD1 HZ]) ■ We take < arg(X) < 2tt, so that the branch 
cut for ^/~X and mi is on the positive real A-axis. 



(6.4) 



The associated spectral density function arising from (|5.9j) is 

MX) :=li ^K(A + ie)] =fc 



1 _ g-Tra/\/A 



Example 2: [q = —a, a < 0; q\ = £(£ + 1)] Repulsive Coulomb 

a ^(£+1) 



y 



y = Xy, a < 0, < x, 00 



2; 

The first Frobenius solution with normalization (|2.9|) is (same as (|6.2|) with a < 0) 

00 

0(x,A) = ^ +1 [1 + £>„(AK] 

n=l 

= x £+1 e ra%/I M(^ + l-ft,2i + 2, -2ixVX) 

= 1 r - jr, M s e+ i (-2ixVX), 

(-2iVX) e+1 +2 



(6.5) 



(6.6) 



where ft := ia/2^/X, and a < 0. The recurrence relation and first three coefficients, a„(A), are the same 
as in Example 1 with a < 0. The Titchmarsh-Weyl m-function arising from (|5.6[) is 



(A) = k t (X) \-alog(-2iVx) - a*(l - ia/(2\A)) - 270 + iVA] + p*(A). 



(6.7) 



where a < 0, and fc^(A), pe(X) are the same, with a < 0, as given above for the hydrogen atom. The 
branch cut is taken again on the positive real A-axis. 



(6.8) 



The associated spectral density function arising from (|5 .9[) is 

ftW = h(X) 



exp(|a|7r/\/A) — 1 



Example 3: [q Q = v 2 - 1/4, v ^ M/2, M = 0, 1, 2 • • • ; qx = Q] Bessel Equation of Non-integer Order 



v 2 - 0.25 



The first Frobenius solution with normalization (|2.6[) is 

{-l) j Xix 2 i 



(x, A) := x 



f+0.5 



2/ = Ay, v>0,v^N/2,N=l,2,.. 



2 v T(v + l)X-»/ 2 x 1/2 J v (VXx) 



{j\{v + l)j2^ 



(6.9) 



(6.10) 
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The Titchmarsh-Weyl m-function arising from (|5.6p is 



7T 



™(A) = , -.v ■ , (6.11) 

2 1 + 1) ■ siniuir) 

where < arg(X) < 2ir, so that the branch cut for \ v is on the positive real A-axis. 
The associated spectral density function arising from (|5.9j) is 

MA) = 22 , +lr A 2 (, + 1) - (6-12) 

Example 4: [go = iV 2 - 1/4, iV = 0, 1, 2, • ■ ■ ;q 1 = 0] Bessel Equation of Integer Order 

-y" + T — y = Xy, a > 0, N = 0, 1, . . . , (6.13) 

a; 

The hrst Frobenius solution with normalization (|2.12p is 



(x, A) := x 



Af+0.5 



2T{N + l)\- N/2 x 1/2 J N {V\x), N = 0, 1,2, • • • . (6.14) 



3- 

The Titchmarsh-Weyl m-function arising from (|5.6j) is 



mo (A) = - log(-2n/A) + 7 - 2^n2 
A w A^-ff 

m ^( A ) = 22j V (jV!) 2 W "( A ) + 2 2A r +l(^! 1 )2 ' N ~ 1 - (6 - 15) 

where < arg(X) < 2n, so that the branch cut for uiq is on the positive real A-axis. 



The associated spectral density function arising from (|5.9|) is 

2 2N + 1 (N\) 

7 Numerical Methods 



M A ) = o2^u»m 2 ' ( 6 ' 16 ) 



In this section and the following two sections we describe some new numerical methods for obtaining 
approximations to the spectral density function (|1.4[) and then the spectral function, by making use of 
the new representation (15.30[) in Theorem 3, and compare their performance with SLEDGE. For general 
information and discussion of numerical methods for Sturm-Liouvillc problems we refer to Pryce's book 
[2Tj and for spectral function computation using SLEDGE we refer to our previous papers [20. M ITS]. 

Many numerical methods for (jl.l[) break down near a singular point at x — 0. However, when we 
take this singular point to be a regular singular point, it admits a convergent Frobenius expansion, and 
then a finite number of leading terms in the sum can be used as an initial approximation near x = 0. 
For equation (|2.4[) the indicial equation is (|2 .5[) and the principal solution is the first Frobenius solution 
with the larger indicial root, 

n = v := 0.5 + 0.5^/1 + 4^0. 

This solution has the general form 

00 

4>{X,X) =J2<*nX n+V , (7.1) 

n=0 

and the general recurrence formula is 

<7ia 
v[y + 1) - q 
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and for n > 1 



— Xa n -2 + q\a n -\ 



n-2 



(1/ + n — l)(z/ + n) — go 

The choice for ao fixes the normalization of (/>(•, A), and in this paper we have made the simple choice 
ao = 1 in all the cases (|2.6|) . (|2.9I) and (12.121) ; this ensures that the properties (i),(ii),(iii) of section 2 
hold in all the cases of Assumption 1. 

We note that there is a risk of loss of significance in the computation of a n (X) for very large A, 
because for moderate n the powers of A in the numerator of a n build up faster than the denominator 
does. We have found that keeping 



x < xo(A) 



\/Vx 



(7.2) 



and using the truncated Frobenius series (to full machine precision) only on (0,Xo(A)] works well. For 
x > xo(A) we use the methods from [12] and [13] for regular problems. In brief, the algorithm is as 
follows: 

(i) For a given A choose a 'matching point' x(A). 

(ii) Use the first Frobenius solution (|7.1I) (this is the principal solution near zero) on the interval (0, xo(A)] 
to produce values for solution and its derivative (usually to machine precision) at xq(A). 

(iii) Apply standard methods to numerically estimate y for (|5.ip on the interval [xo(A),x(A)] satisfying 
initial conditions from (i). Because of the oscillation in y for A > 0, we use a piecewise trigonometric 
approximation to y. 

(iv) Approximate the solution P(x, A), Q(x, A), and R(x, A) of the initial value problem (|1.6p - (|1.7[) at 
the matching point x(A) using one of the approaches described below. 

(v) Substitute the estimates from (iii) ,(iv) into 



(7.3) 



tt[P(x, A)y(x, A) 2 + Q(x, X)y(x, X)y'(x, A) + R(x, \)y'(x, A) 2 ] ' 

to produce an approximation to the exact spectral density, /(A), given in Theorem 3, equation (|5.30l) . 

The two papers [12] and [13] derive two very different approaches to the numerical computation of 
(P(b, A), Q(b,X), and R(b, A)). In [12] we constructed a family of recurrence formulas that generated 
successively more accurate approximations to P, Q, and R and hence to /(A). For fixed x > define 
for each positive integer j the family of functions 



nlP^ + Qjyy' + Rjy' 2 ]' 



From [TT] for j = 1 we define 



Pi := VA, Qi := 0, Pi := 1/V\. 
From [12] the next formula in the family for j = 2 is defined by 

P 2 := y/\-q(x), Q 2 := -?'(x)/[2(A - q(x))^% R 2 := 1/y/X - q(x). 
The final member that we use is j = 3, also defined in 



as: 



P 3 

Q 3 
R 3 



P 2 + 0.25 7 2 + 0.125 7l 2 /7o 
P 2 - 0.25 7 2 + 0.125 70 7 2 



where 



Ik := 



dx k 



(7.4) 



(7.5) 



(7.6) 



(7.7) 
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for k = 0, 1, 2. For regular problems on [A, oo), A > 0, it is shown in [12] that each member of this family 
converges to the spectral function /a(A) as x — > oo. For the hydrogen atom potential (|6.1|) on [.A, oo) 
the theory of 12 implies that 



/.4(A) - F*(A) = Oil/x 2 *- 1 ) as z -> oo 



(7.8) 



This method requires knowledge of derivatives of the potential q(x). 

In [13] we constructed explicit approximations to the solutions of (|1.6|) with known residual terms 
that arise during the construction; in particular, replace P(x, A), Q(x, A), and J2(x, A) with estimates of 
the form 

N 

P N (x) := VA + ^2 a iM 

3=1 

TV 

Qjv(s) := E 6 ^' +1 
3=1 

N 

R N (x) := l/Vx + ^Cj/x*, 

3=1 



(7.9) 



where {fly}, {^j}, and {a,} will depend on A but not x. The resulting sums are substituted into (|1.6[) 
and the coefficients chosen to match terms of like powers. Then we put the computed coefficients into 
(|7.9|) and define the family of approximations 



/^(A) 



1 



n[P N (x, X)y(x, A) 2 + Q N (x, X)y(x, \)y'{x, A) + R N (x, X)y'(x, A) 2 ] 
Specifically, the iVth residuals are defined as 



(7.10) 




P 1 



Q'i 

R'r 



N 



N 



X-q 
-2 2(A-q) 
0-1 



Pn 
Qn 
Rn 



and we attempt to make these small, as x — > oo, by the choice of coefficients in (|7.9[) . All potentials in 
the examples of the previous section have the form ()5.5|) . that is, 

(7.11) 



q(x) = A/x + B/x 2 , 

where A = q\ and B = qo and (|2.2[) is satisfied. It is straightforward to show that 

N 



N 



-J a 3 



Xbj + Abj-x + Bbjs 



r-3 + l 



Bb N -i + Ab 



„N+2 



E 

3=1 
N 

E 

3=1 

Nb N _ 1 + 2Ac N + 2BCN-X 2Bc N -(N + l)b 



N 



Bb 



N 



JV+3 



-{j - l)bj- 2 + 2aj - 2Xcj + 2Acj-! + 2Bc d - 2 2A + 2B/x 



cVx 



N 



„N+2 



4 



N 

E 

3=1 



~3+l 



If we require the coefficients to satisfy 

jdj + Xbj = Abj-i + Bbj-2 

= (j - l)6j_ 2 /2 - Acj-i - Bcj-2 + [ASp + B6 j2 ]/Vx 



2j - XCj 



jcj = 0, 



(7.12) 
(7.13) 
(7.14) 
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for j — 1, 2, . . . , N, then the residuals simplify to 



P _ Ab N + Bb N -i Bb N 

»N - h x N+3 V- Lt >) 

,q —Nbfif—i + 2Acn + 2Bcn-i 2Bc N -{N+l)b N 

b N = -pm + ( 7 - 16 ) 



0. (7.17) 



If we adopt the convention that coefficients with nonpositive subscripts have zero values, then the solution 
of ([7TT2j) - (fTT3]l can be written, for 1 < j < N, 

aj = (h+t 2 )/2 
Cj = (t 1 -t 2 )/(2X) 
bj = jCj, 

where 

t x = [(Abj^ + Bbj^/j 

and _ 
*a - 0.5(j - l)6j_ 2 - Acj-x - Bcj-2 - [A6 n + B5 j2 ]/Vx. 

Since the derivatives of the residuals do not change sign once x is sufficiently large, the theory of [T!2] 
implies for that 

f(\)-f?(\)=0(l/x N+1 ) (7.18) 

as x — > oo. 

To numerically estimate the spectral density function /(A), we would usually use the methods (|7.4[) 
from [TJ] because they require knowledge of only the first few derivatives of q(x). But when q has the 
required special forms, the method (|7.10p from [13], often more efficient, can also be used. 

8 Numerical Estimation of the Spectral Density Function /(A) 

In this section we test our implementation of the various numerical methods from the previous section 
on the examples listed in section 6, for which exact formulas are known for the spectral density function. 
Then we also test a more interesting example from quantum chemistry for which exact formulas are not 
known. 

For the hydrogen atom potential, Example 1 (equation (|6.ip ) with a = 1, Table 8.1 has numerical 
output when I = 1 for many A with the methods F^, F%, F%, and f% using the notation of the previous 
section. Table 8.2 shows the analogous data when 1 = 2. Shown are only the errors: absolute when 
the answer is less than one and relative otherwise. A tolerance of 10 -14 was used for the numerical 
integration of the initial value problem for (|6.1[) . starting at ieo(A) from (|7.2|) . Consequently, table 
entries this small represent errors in y as well as errors due to finite x. Note that for a fixed accuracy, 
generally larger matching points x are needed when A is smaller. As expected, the higher order methods 
F 3 and / 6 are much superior. 



25 



Table 8.1. Numerical Error: H Atom (£ = 1) q(x) = -l/x + 2/x 2 . 



A 


x = x(a) 


F 




F 




F 




fx 




0.1 


320.0 


4.60( 


-4) 


-2.10( 


-8) 


5.62( 


-12) 


2.13(- 


14) 


0.2 


225.0 


2.19( 


-5) 


-1.18( 


-9) 


4.28( 


-13) 


4.48(- 


14) 


0.4 


160.0 


5.08( 


-4) 


— 2.34( 


-8) 


6.36( 


-12) 


4.56(- 


15) 


1 


100.0 


5.53( 


-4) 


— 2.62( 


-8) 


7.49( 


-12) 


1.25(- 


13) 


2 


71.0 


9.52( 


-4) 


-4.39( 


-8) 


1.21( 


-11) 


9.60(- 


14) 


4 


50.0 


-2.55( 


-4) 


1.22( 


-8) 


-2.95( 


-12) 


4.89(- 


13) 


10 


32.0 


1.41( 


-3) 


-5.93( 


-8) 


1.48( 


-11) 


3.25(- 


13) 


20 


22.5 


-2.79( 


-4) 


1.14( 


-8) 


-1.81( 


-12) 


7.83(- 


13) 


40 


16.0 


3.52( 


-4) 


-1.21( 


-8) 


2.91( 


-12) 


8.16(- 


13) 


100 


10.0 


-3.39( 


-4) 


8.44( 


-9) 


2.52( 


-13) 


2.61(- 


13) 


200 


7.0 


2.32( 


-4) 


-2.34( 


-9) 


-1.91( 


-12) 


3.02(- 


13) 


400 


5.0 


-1.04( 


-4) 


-1.66( 


-9) 


2.78( 


-121 


2.60(- 




1000 


3.2 


-4.06( 


-6) 


-6.53( 


-10) 


7.80( 


-13) 


2.82(- 


13) 


2000 


2.2 


5.84( 


-6) 


5.12( 


-9) 


-2.91( 


-12) 


2.86(- 


13) 


4000 


1.6 


3.49( 


-6) 


-1.79( 


-9) 


1.34( 


-12) 


2.84(- 


13) 


10000 


1.0 


2.67( 


-5) 


-6.60( 


-9) 


3.81( 


-12) 


2.95(- 


13) 


Table 8.2. Numerical Error: 


H Atom (£ = 




—1/x 


+ 6/x 2 




A 


x = x(X) 


F 1 

* X 




F' 2 




F 

X 




fx 




0.1 


320.0 


-3.13(- 


-6) 


1.35(- 


-10) 


-3.45(- 


14) 


3.95( 


-15) 


0.2 


225.0 


-4.79(- 


-6) 


2.13(- 


-10) 


-5.60(- 


14) 


2.98( 


-15) 


0.4 


160.0 


-4.28(- 


-6) 


1.83(- 


-10) 


-4.54(- 


14) 


3.00( 


-15) 


1 


100.0 


-2.82(- 


-5) 


1.20(- 


-9) 


-2.98(- 


13) 


8.32( 


-15) 


2 


71.0 


-7.24(- 


-5) 


2.87(- 


-9) 


— 6.41(- 


13) 


1.92( 


-14) 


4 


50.0 


8.09(- 


-5) 


-2.97(- 


-9) 


5.83(- 


13) 


5.28( 


-15) 


10 


32.0 


-8.59(- 


-4) 


2.24(- 


-8) 


-7.48(- 


13) 


1.69( 


-13) 


20 


22.5 


2.80(- 


-4) 


-3.67(- 


-9) 


— 1. / o\ — 


1 9\ 
11) 




-16) 


40 


16.0 


-2.36(- 


-4) 


-2.06(- 


-9) 


4.70(- 


12) 


6.52( 


-14) 


100 


10.0 


1.69(- 


-4) 


1.66(- 


-8) 


-1.22(- 


11) 


8.75( 


-14) 


200 


7.0 


-4.57(- 


-5) 


-2.62(- 


-8) 


1.71(- 


11) 


3.02( 


-14) 


400 


5.0 


-3.42(- 


-5) 


2.25(- 


-8) 


-1.26(- 


11) 


-5.91( 


-14) 


1000 


3.2 


-8.80(- 


-6) 


2.46(- 


-9) 


-6.75(- 


13) 


8.21( 


-13) 


2000 


2.2 


1.06(- 


-4) 


-2.25(- 


-8) 


1.24(- 


11) 


8.62( 


-13) 


4000 


1.6 


-2.78(- 


-5) 


5.40(- 


-9) 


-2.15(- 


12) 


8.09( 


-13) 


10000 


1.0 


-1.42(- 


-4) 


2.38(- 


-8) 


-1.17(- 


11) 


8.94( 


-13) 



The next choice is the Bessel Equation, Example 3 (equation (|6.9p l. with v = 1/3: 

q(x) = -5/(36a: 2 ). 
The error behavior is similar to that in the previous tables. 
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Table 8.3. Numerical Error: q(x) = -5/(36a; 2 ) 



A 


/ \ \ 

x = x(a) 


r 

X X 




7^2 

b 

^ X 




r 

^ X 




fx 




0.1 


320.0 


— 3.26( 


-7) 


4.61(- 


-11) 


3.21(- 


-14) 


5.37(- 


14) 


0.2 


225.0 


1.22( 


-6) 


-1.13(- 


-10) 


1.59(- 


-13) 


6.86(- 


14) 


0.4 


160.0 


-5.18( 




7.31(- 


-11) 


4.90(- 


-14) 


8.34(- 


14) 


1 


100.0 


2.34( 


-6) 


— 3.48(- 


-10) 


2.85(- 


-13) 


1.12(- 


13) 


2 


71.0 


3.30( 


-6) 


-4.93(- 


-10) 


3.92(- 


-13) 


1.47(- 


13) 


4 


50.0 


3.71( 


-6) 


-5.53(- 


-10) 


4.38(- 


-13) 


1.63(- 


13) 


10 


32.0 


-1.51( 


-6) 


2.14(- 


-10) 


1.49(- 


-13) 


2.49(- 


13) 


on 


99 z. 




(\\ 
-v) 




-10) 


0. ( o^ — 


-16) 


9 8^^ 


16) 


40 


16.0 


-1.78( 


-6) 


2.51(- 


-10) 


1.58(- 


-13) 


2.76(- 


13) 


100 


10.0 


5.92( 


-6) 


-8.26(- 


-10) 


7.22(- 


-13) 


2.83(- 


13) 


200 


7.0 


-5.92( 


-6) 


9.12(- 


-10) 


-1.74(- 


-13) 


2.95(- 


13) 


400 


5.0 


5.92( 


-6) 


-8.82(- 


-10) 


7.36(- 


■13) 


2.97(- 


13) 


1000 


3.2 


-1.78( 


-6) 


2.51(- 


-10) 


1.54(- 


-13) 


2.72(- 


13) 


2000 


2.2 


-5.78( 


-6) 


8.89(- 


-10) 


-1.65(- 


■13) 


2.91(- 


13) 


4000 


1.6 


-1.78( 


-6) 


2.51(- 


-10) 


1.70(- 


-13) 


2.88(- 


13) 


10000 


1.0 


5.92( 


-6) 


-8.82(- 


-10) 


5.98(- 


-12) 


2.63(- 


13) 



In the papers of Bain et al [4] , Brandas et al [5] , and Engdahl et al [7] , [8] , [9] can be found a potential 
with a "barrier" near x = 2 and decaying rapidly to zero as x — > oo: 



. , £(£ + 1) a ir 2 _ x 
q(x) = 9 h 15x 2 e x . 

x x 

For this example we have qo = 1(1 + 1), q\ = —a, qi = qj, = 0, and for k > 2 

15 



(8.2) 



9fe+2 



.if 



(fc-2)f 



For a = 1, Table 8.4 displays the results for several values of x to show the rapid convergence as x — > oo. 
A tolerance of 10 -14 was used for the numerical integration of y. 



Table 8.4. Estimates for Barrier Potential (IQl with a = 1 



X 




£ = Q 




1=1 




1=2 




e = o 


e = i 


£ = 2 










A = 7 










A = 10 




5. 





142809355 





019804657 





0004166628 


1 


686525464 


1.728228916 


0.1242771047 


10. 





142828980 





019801387 





0004162081 


1 


686646374 


1.728086680 


0.1242724151 


15. 





142829143 





019801395 





0004162075 


1 


686647533 


1.728085796 


0.1242722883 


20. 





142829149 





019801396 





0004162075 


1 


686647559 


1.728085772 


0.1242723271 


25. 





142829149 





019801396 





0004162075 


1 


686647559 


1.728085772 


0.1242723271 










A = 20 










A = 40 




5. 


1 


999374819 


4 


314112367 


2 


8784054457 


2 


558971562 


11.31991454 


17.270757428 


10. 


1 


999374881 


4 


314112204 


2 


8784044730 


2 


558971266 


11.31991563 


17.270756351 


15. 


1 


999374882 


4 


314112311 


2 


8784043230 


2 


558971293 


11.31991552 


17.270756523 


20. 


1 


999374882 


4 


314112307 


2 


8784043239 


2 


558971293 


11.31991552 


17.270756528 


25. 


1 


999374882 


4 


314112307 


2 


8784043240 


2 


558971293 


11.31991552 


17.270756528 



9 Numerical Estimation of the Spectral Function p(X) 

Associated with the density function /(A) is the associated spectral function defined by 

p(X) := C f(p)dii. 



(9.1) 
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In [TT] for problems regular at x — we estimated p using F 1 and compared with the package SLEDGE 
[201 US] ; generally the p(X) computation using a quadrature routine for (|9.1j) and the F 1 formula ran 
considerably faster than SLEDGE, but still had the drawback that rather large a;— intervals were required 
for the F 1 calculation. Here we apply the methods of this paper for computing p by estimating / in 
(|9.1[) . and performing a quadrature, to the examples in section 6 for which exact answers are known, 
and again compare with SLEDGE. 

The SLEDGE software for estimating p(X) is based on the Levitan-Levinson characterization of the 
spectral function as a limit of step spectral functions over a finite interval approximation (second formula 
in (15.81) ); this is a totally different approach than the present approach of this paper which relics on the 
family of i^-approximants, together with the quadrature in (|9.1[) . For the case of two singular endpoints, 
the performance of SLEDGE for computing the spectral function on examples having explicit formulas 
for the spectral function was reported on in [15]. As reported there, one of the major weaknesses of 
the SLEDGE package is obtaining high accuracy in the p(X) calculation when A is large; this is due 
primarily to the fact that SLEDGE does not rely on asymptotic approximations for the eigenvalues and 
eigenfunction norm reciprocals, but computes them numerically as required for implementing the Pb(X)- 
formula in (|5.8I) . Experience in using SLEDGE on doubly singular problems is that very large computing 
times are required due to the computation of large numbers of eigenvalue - eigenfunction norm pairs, 
and that there is significant loss of accuracy when A becomes sufficiently large. As the timing and 
accuracy data of this section shows, doubly singular problems can be handled with high accuracy and 
much reduced computing times by making use of the F^-approximants and the /^-approximants of this 
paper, along with the quadratures for computing p(X) in (|9.ip ; this represents a major improvement in 
computational technique over the SLEDGE algorithm for spectral function computation. 

Following SLEDGE, we assume approximations are sought for a finite set of A- values in the continuous 
spectrum, (0, oo), ordered so that 

< Ai < A 2 < ... < A m . 
Then with p(0) given (or computed via SLEDGE) we estimate 

p(X 1 ) = p(0)+ [ 1 f( f i)dp 
Jo 

and for j = 2, 3, . . . , m 

p(A i )=p(A i _ 1 )+ / 3 f(n)dn (9.2) 

using an adaptive quadrature code with / replaced by FJ or approximations. Here we report some 
timing and accuracy data for the four examples listed in section 6. 

The spectral functions on (0, oo) for these examples are known in closed form by putting the exact 
spectral density functions from (|6.4p . (|6.8I) . (|6. 12|) . and (|6.16|) into (|9 . 1 1) and performing an exact inte- 
gration. The resulting closed form formulas for p(X) were used for the four examples to compare with 
the numerical approximations and to generate the 'exact' error; the errors are taken as absolute if the 
exact value of p(X) is less than one, and relative otherwise. 

For the Bessel equation of order 1, Example 4 (equation (16.131 with N=l), we used N = 7 in the 
approximation (|7. 10[) (that is, the scheme from [T3l Sec 4]). Output data for six A- values is displayed in 
Table 9.1. The quadrature tolerance was 10 -8 and the tolerance for the initial value problem was 10 -9 . 
Note that, at this tolerance, all the apparent error arises from the first integration interval and is passed 
on through the sum in (|9.2|) . 



Table 9.1. p(A) for the first order Bessel equation on (0, oo) (q(x) = 0.75/x 2 ) 



X 


A = 1 


A = 2 


A = 4 


A = 10 


A = 20 


A = 40 


6. 


0.06318042 


0.25068042 


1.00068041 


6.25068041 


25.00068041 


100.00068041 


12. 


0.06254252 


0.25004252 


1.00004252 


6.25005252 


25.00004252 


100.00004252 


24. 


0.06250266 


0.25000266 


1.00000266 


6.25000266 


25.00000266 


100.00000266 


36. 


0.06250000 


0.25000000 


1.00000000 


6.25000000 


25.00000000 


100.00000000 
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The data for estimating /(A) in section 8 showed that as A gets larger, smaller values of the matching 
point x — x(X) are needed for a given accuracy. The data in Table 9.1 for p exhibit this phenomena. 
Hence, from efficiency considerations in order to compute p(X) we want the choice of matching point x 
to vary with the integer N, or even better, pointwise with A. For q(x) — 0.75/ x 2 it can be shown that 
the absolute error in /(A) for a given N is proportional to 

1 

^+2^-1/2- 

This suggests that taking x ~ \/ \( 2N - 1 )/( 2N + 4 ) would be a good heuristic for the matching point. 
Similarly, for the family of approximants, F , from 1|7.4|) — (|7.7p the corresponding form for the absolute 
error in /(A) is 

1 

x 2N X N - 1 ^' 

so that x ~ \ / \( 2N - 1 ) / ( iN ) would be an appropriate matching point. The latter is roug hly 1/y/X. As 
mentioned earlier, these formulas would change for a different q. Also, the heuristic would differ for 
relative errors. A similar analysis for the general potential q(x) = A/x + B/x 2 suggests a good first 
value of matching point would be 

x = x{\) := \A\/(2X) + v/A 2 /(4A 2 ) + \B\/\. (9.3) 

We have written a research code, called AutoB, for which the only inputs required are the set of 
A points, the p(0) value, and the accuracy desired. If q has the form of (|7.11l) or the form of similar 
potentials in [TJ], then we use the appropriate formula in (|7.10p with N chosen to be a function of 
the accuracy sought. Otherwise, we use F 3 from (|7.4I) which requires knowledge of derivatives of q. We 
report the performance of AutoB on the four examples in section 6 using (|9.3[) as the initial choice of 
matching point. Given a prescribed tolerance r, for each A the matching point x is then increased until 

(error estimate| < max{l, |output value|} r 

holds. Estimates at various r were sought for the following set of sixteen A values: 

{0.1, 0.2, 0.4, 1, 2, 4, 10, 20, 40, 100, 200, 400, 1000, 2000, 4000, 10000}. (9.4) 

The error shown is the maximum (relative when the 'exact' p > 1, absolute otherwise) over the set 
of sixteen A values. For many of these runs the heuristics were overly conservative, but the times are 
nevertheless quite small. For the Bessel equation of order ^, Example 3 (equation (|6.9[) with v — 1/3), 
the relatively large computing times were due to difficulties near A = 0. 

Table 9.2. AutoB results for several tolerances and Four Potentials on (0,oo). 







T = 10" 


-4 


T = 10" 


-6 


T = 10- 


-8 


T 


= io- 


-10 


Potential 




error 


time 


error 


time 


error 


time 


error 


time 


1. Ex4(i/= 


=1) 


1.30(-7) 


0.28 


1.50(-9) 


0.50 


1.92(-12) 


1.06 


1.73(- 


-13) 


2.46 


2. Exl(£= 


=1) 


1.78(-6) 


0.42 


7.88(-7) 


0.84 


1.50(-8) 


3.11 


7.93(- 


-11) 


18.50 


3. Ex2(£= 


=1) 


2.49(-6) 


0.37 


3.37(-7) 


0.74 


8.54(-9) 


3.95 


1.10(- 


-10) 


20.97 


4. Ex3(^= 


=1/3) 


9.57(-5) 


0.49 


1.04(-6) 


3.00 


1.52(-8) 


41.95 


1.13(- 


-10) 


894.77 



For comparison, output from the SLEDGE program is shown for the Bessel equation of order Example 
4 with v = 1, Example 4 (equation (|6.13[) in Table 9.3. Since SLEDGE is known [T5] to have difficulty 
with large values of A, we ran the program at various choices of r only on the first n A- values from 
(|9.4j) with n — 1, 2, ... ,9. The final line (n = 16) is output from AutoB using all sixteen A values up 
to A = 10000. Clearly, AutoB is much more reliable than SLEDGE. Similar results were observed for 
other doubly singular potentials. 
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Table 9.3. SLEDGE output for the Bessel equation of order l(q(x) = 0.75/x 2 ) 







T = 10 


— 3 


T = 10 


-4 




r = 


10 -5 


r = 


10 -6 


n 


error 


time 


error 


time 


error 


time 


error 


time 


1 


1.31 


-4) 


0.02 


3.87(— 5) 


0.09 


9.74( 


-6) 


0.50 


9.92(— 6) 


1.58 


2 


6.86 


-4) 


0.02 


1.24(-5) 


0.23 


1.07( 


-5) 


1.31 


9.87(-6) 


8.43 


3 


2.72 


-4) 


0.14 


1.32(-5) 


0.34 


2.95( 


-5) 


2.34 


6.70(-6) 


>96.53 


4 


2.71 


-4) 


0.16 


1.79(-4) 


0.50 


3.20( 


-5) 


8.45 


6.70(-6) 


>177.28 


5 


7.16 


-4) 


0.19 


1.79(-4) 


0.70 


3.20( 


-5) 


16.40 






G 


1.92 


-3) 


0.21 


1.79(-4) 


1.71 


3.20( 


-5) 


29.51 






7 


7.04 


-3) 


0.96 


1.83(-4) 


7.67 


8.65( 


-5) 


51.64 






8 


1.47 


-2) 


1.83 


1.83(-4) 


11.39 


2.18( 


-4) 


75.87 






9 


2.23 


-2) 


>12.00 


3.31(-4) 


20.48 


2.83( 


-4) 


>465.99 







AutoB 



16 3.69(-6) 0.18 1.30(-7) 0.28 5.74(-9) 0.36 1.50(-9) 0.50 

To illustrate the superiority of the new code AutoB over SLEDGE we also ran comparisons on timing 
and accuracy the Hydrogen Atom potential with I = 1, Example 1 (equation (|6.1j) . The output values 
for p(A) obtained for each of the sixteen A values in (|9.4|) for four choices of the tolerance levels are 
displayed in Table 9.4. The > in the time needed for SLEDGE indicates that it stopped (too much 
time) before the user requested input accuracy was achieved. Since SLEDGE could not achieve 10 -3 - 
accuracy over the whole range of A-values, no SLEDGE runs for tighter tolerances are listed. As the 
data shows, SLEDGE has much difficulty to compute highly accurate results for large values of A, while 
the new codes are capable of quite high accuracy in much less computing time. Similar testing for the 
Hydrogen Atom potential using the approximants was also done in the thesis of Mark Schuster [22] . 



Table 9.4: Comparison of SLEDGE with AutoB for Hydrogen problem with 1=1. 



A 


Exact 


SLEDGE 


AutoB 


AutoB 


AutoB 


AutoB 


0.1 


0.005621362 


0.0056 


0.0056 


0.00562 


0.0056214 


0.005621362 


0.2 


0.010067470 


0.0100 


0.0100 


0.01007 


0.0100675 


0.010067470 


0.4 


0.022334469 


0.0222 


0.0223 


0.02233 


0.0223345 


0.022334470 


1.0 


0.087358065 


0.0867 


0.0874 


0.08736 


0.0873581 


0.087358074 


2.0 


0.298717032 


0.2966 


0.2987 


0.29872 


0.2987170 


0.298717047 


4.0 


1.166166722 


1.1577 


1.1662 


1.16617 


1.1661667 


1.166166736 


10. 


8.206942681 


8.1493 


8.2069 


8.20694 


8.2069402 


8.206942643 


20. 


38.98117554 


38.691 


38.981 


38.9812 


38.981169 


38.98117542 


40. 


194.5884791 


192.80 


194.59 


194.588 


194.58845 


194.5884785 


100. 


1719.215348 


1706.1 


1719.2 


1719.21 


1719.2140 


1719.215343 


200. 


9188.295022 


9068.2 


9188.3 


9188.29 


9188.2922 


9188.294986 


400. 


49923.13741 


49137. 


49923. 


49923.1 


49923.126 


49923.13720 


1000. 


475962.2250 


462178. 


475962. 


475961. 


475962.23 


475962.2484 


2000. 


2644112.132 


2479326. 


2644120. 


2644111. 


2644111.6 


2644112.120 


4000. 


14766762.30 


13806576. 


14766758. 


14766759. 


14766759. 


14766762.23 


10000. 


144274264.9 


112742503 


144274122 


144274122 


144274123 


144274264.3 


time (sec) 




>177.5 


0.31 


0.42 


0.83 


3.75 


RelErr 




10^3 


10-3 


io- 4 


io- 6 


io- 8 


AbsErr 




10^3 


10-3 


io- 4 


io- 6 


io- 8 



Remark. High accuracy in the spectral function computation for the Bessel equation on (0,oo), 
Examples 3 and 4, was also achieved in [TSl Sec 4]; this, however, was done by inserting asymptotic 
formulas for the eigenvalues and eigenfunction norm reciprocals for the Bessel equation on (0, b] into the 
SLEDGE code (bypassing the SLEDGE computation of these quantities); but, of course, this was not 
an automatic procedure applicable to other problems with two singular cndpoints. 
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